Fragmentation and Evolution of Molecular Clouds. Ill: The Effect of Dust and 

Gas Energetics 

Hugo Martel 1 ' 2 , Andrea Urban 3 , and Neal J. Evans II 4 
ABSTRACT 

Dust and gas energetics are incorporated into a cluster-scale simulation of star for- 
mation in order to study the effect of heating and cooling on the star formation process. 
We build on our previous work by calculating separately the dust and gas temperatures. 
The dust temperature is set by radiative equilibrium between heating by embedded stars 
and radiation from dust. The gas temperature is determined using an energy-rate bal- 
ance algorithm which includes molecular cooling, dust-gas collisional energy transfer, 
and cosmic-ray ionization. The fragmentation proceeds roughly similarly to simulations 
in which the gas temperature is set to the dust temperature, but there are differences. 
The structure of regions around sink particles have properties similar to those of Class 
objects, but the infall speeds and mass accretion rates were, on average, higher than 
those seen for regions forming only low-mass stars. The gas and dust temperature have 
complex distributions not well modeled by approximations that ignore the detailed ther- 
mal physics. There is no simple relationship between density and kinetic temperature. 
In particular, high density regions have a large range of temperatures, determined by 
their location relative to heating sources. The total luminosity underestimates the star 
formation rate at these early stages, before ionizing sources are included, by an order 
of magnitude. As predicted in our previous work, a larger number of intermediate mass 
objects form when improved thermal physics is included, but the resulting IMF still 
has too few low mass stars. However, if we consider recent evidence on core-to-star 
efficiencies, the match to the IMF is improved. 

Subject headings: hydrodynamics — ISM: clouds — ISM: dust — methods: numerical 
— stars: formation 

1. Introduction 

The physics of star formation is the link between the small-scale - planet formation - and the 
large-scale - galactic evolution. Many, probably most, stars form in highly clustered environments 
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(|Lada fc Ladal 120031 ; iBressert et al.l l2010l ) . Within the highly complex stru cture of a molecular 
cloud, theorists have ide ntified the clump as the object that forms a cluster ([Williams et al.ll2000i ; 
McKee &: Ostrikerl 120071 ) . Observational studies have found a range of structures that might cor- 
respond to this concept, but the formation of rich clusters seems most clearly associ ated with 



parti cularly dense clumps, identified by strong emission from tracers of dense gas (e.g., IWu et al 
2O10h . 



This paper is the third of a series that studies the fragmentation of a dense molecular clump 
using Smoothed Particle Hydrodynamics (SPH ) to follow the hy drodynamics, with a focus on 
the effects of the thermal physics. In Paper I (jMartel et al.ll2006l . hereafter MES06), we showed 
that an isothermal gas will fragment excessively, producing only very low mass stars. In Paper II 
(jUrban et al.ll2010l . hereafter UME10), we included global radiative feedback from the forming stars, 
assuming that the gas temperature was equal to the dust temperature; this approximation over- 
produced high mass stars. In this paper, we calculate the gas temperature separately from the dust 
temperature. We use these simulations to address the role of thermal physics in the fragmentation 
problem, the density distribution around forming stars in a proto-cluster, the evolution of the far- 
infrared luminosity during formation, with application to the use of far-infrared luminosity as a 
probe of star formation rate, and the effects of an improved thermal physics on the mass function 
of forming stars. 

In our previous work (UME10) we modeled a clustered star-forming region and showed that 
dust-gas thermal energetics with source luminosity terms from young stars can heat the gas and 
prevent fragmentation. Less fragmentation led to the formation of massive stars, which had not 
been produced in our isothermal simulation (MES06). Although we were able to form massive stars, 
we missed a significant fraction of the low-mass stellar population. We hypothesized in UME10 that 
by including a more realistic dust-gas thermal energetics algorithm we would increase the number 
of low-mass stars. Including molecular cooling would decrease the temperature of the gas, leading 
to more fragmentation and more low-mass stars. To t est our hypothesis, we have implemented 
the complete heating and cooling algorithm described in I Urban et al.l ( 2009 ) (hereafter UED09) in 
simulations with similar scales and parameters as the simulations discussed in UME10. We discuss 
our work in the following sections. In $2j we discuss our numerical algorithm and our new method 
of calculating the gas temperature. In we discuss the initial conditions and parameters of our 
simulation. Our results are discussed in £JH We conclude and summarize the paper in £}5j 



2. The Numerical Algorithm 



Our numerical algorithm was described in MES06 and U ME10. It is a standard Smoothed 
Particle Hydrodynamics (SPH) algorithm (see lMonaghanlll992l . references therein), which simulates 
the growth of structures in a cubic volume with periodic boundary conditions, representing a 
small part of a giant molecular cloud. The code was modified to include particle splitting and 
sink particles. In the optically-thin regime, the Jeans mass decreases with increasing density. 
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Eventually, when the density becomes sufficiently high, the gas becomes optically thick and the 
Jeans mass starts increasing with density. Hence, there is a minimum Jeans mass, corresponding 
to the transition between the two regimes. To properly follow the fr agmentation of the cloud, it is 



essential to resolve that minimum Jeans mass. Particle splitting (see iKitsionas &: Whitworthll2002l ; 
MES06) enables us to do this at a reasonable computational cost. 

When the gas reaches a c ertain critical density p c , the algorithm replaces gas fragments by sink 



particles, using the method of lBromm et al.l (|2002h . Sink particles (or sinks) represent protostellar 



cores. They are not allowed to fragment or merge, but they have the ability to grow by accreting 
surrounding gas particles. Any bound gas particle within its accretion radius, r acc (~ 150 AU 
for all but one simulation, see Table [1]), is automatically accreted into the sink. Because r acc is 
considerably larger than the actual forming star (likely radius a few to tens of solar radii), the 
evolution of material inside the sink is unknown. For the simulations, we assume that all the 
material falling into the sink flows continuously onto the actual stellar core, producing accretion 
luminosity. The luminosity resulting from that accretion will heat the surrounding gas. In order t o 



calculate the luminosity from a sink particle, we use the models of lWuchterl &: Tscharnuterl (|2003l ). 
specifically their Table 3, which include the effect of mass accretion on the luminosity. For objects 
with masses greater than 2M Q , we use the method described in UME10 to calculate the luminosity. 
There is considerable evidence that not all material passing through a radius of 150 AU winds up 
in the forming star, so we also consider the effects of the loss of some material in comparing to the 
mass function fiJ4.5D. 



2.1. Dust and Gas Temperature Calculation 

We use the same numerical methods described in UME10 to calculate the mass accretion rate 
onto sinks and the dust temperature. In order to calculate the gas temperature, we use the method 
of UED09. We give a brief description here. 

The dust temperature is determined using the method discussed in UME10, an d in more detail 



in UE D09. UED09 used a spherically-symmetric radiative transfer code, DUSTY (jNenkova et al 



2000), to calculate the dust temperature distribution around young stellar objects. Using DUSTY, 
we created a grid of models with input values of luminosity and density distribution. In the 
simulations presented in this paper, we calc ulate the luminosity of individu al objects based on mass 



and mass accretion rate using the models of lWuchterl fc Tscharnuterl (|2003l ) . We also determine the 



density profile around each of the sinks formed in our simulation. We then use the grid calculated 
in UED09 to find an analytic fit to the dust temperature distribution around each of the individual 
sinks. We approximate the dust/gas distribution around sink particles as spherical. 



We calculate the density profile around individual sink particles, as in UME10, using spherical 
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shells. We parameterize the density profile with n Q and a, as the following, 



n(r) = n Q ^ 



1000 AU 



cm 



(1) 



The symbol n represents the number density of all particles (n = 
ratio nn 2 /nHe = 5, which corresponds to a mean molecular weight pL 
p = fj,m K n. 



n H2 + n 



He) 



We assume a 



2.33. The gas density is 



Our approximations are complementary to those of most other simulations t hat include some 
thermal physics. We do not include compressional heating during collapse, as does lBatd (|2009i ). We 
are no t doing radia t ive tr ansfer during the SPH calculation, but instead using a pre-computed grid 
as did ISmith et al.l (|2009i ). and we are assuming spherical distributions of material around sinks. 
Thus, we will miss some of the effects of geometry included i n papers that do not assume spherical 
symmetry (e.g. , iKrumholz et al-lbooj krumholz et alJboiol . lBatj2009l . loffner et ailhood ) . On the 
other hand, we used realist ic grain opacitie s as a function of wav elength (i.e., OH5 dust opacities, 
Ossenkopf &: Henningll 19941 . as described in lYoung &: Evansll2005l ) in radiative transfer calculations 
that include non-isotropic scattering and apply to all relevant optical dept hs, in contrast to mean 



opacities and other a pproximations used by most other simulations (e.g., IKrumholz et al. 
Krumholz et al.lboij v 



2007 



The gas temperature algorithm (UDE09) includes energy transfer between gas and dust via 
collisions, gas heating by cosmic rays, and molecular cooling. Heating by photoelectric emission 
from dust grains is not included because the external ultraviolet radiation is strongly attenuated 
and the clump we are simulating is assumed to reside deep in a larger molecular cloud. We are 
not includin g ultraviolet radiation from the forming stars. The dust to gas ratio is taken to be 



4.86 x 10~ 3 (IHollenbach fc McKeelll989l ). and the grain cross section p er baryon is 6.09 x 10 
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( Young et al.ll2004r ). The cosmic ray ionization rate is 3 .0 x 10 17 s 1 (van der Tak &: van Dishoeck 



2001 



The fractional abundance 



2000) and the energy deposited per ionization is 20 eV (jGoldsmithl l! 
of CO is taken to be 10 -4 . The algorithm requires inputs of dust temperature (discussed in 
the previous paragraph), local velocity dispersion (5v), column density, and local density. We 
use the local density calculated from the density profile [eq. ([TJ] to be self-consistent with our 
dust temperature calculation. The local velocity and the column density are needed in order to 
calculate the level of radiative trapping in the molecular cooling lines. The local velocity dispersion 
(characterized by the Doppler b parameter) was assumed to be 1 km s _1 throughout the calculation. 
This is a reasonable assumption based on the values of the velocity field found in UME10 (see Table[5] 
below). 

We estimate the column density at a point of interest r within our simulation with the following 
line integral, for every sink i present in the simulation at that time: 



N, 



column, i 



A«i+2000AU 



Al, 



rii(r)dr, 



(2) 
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where AZj is the distance from sink % to r, nj(r) is the density profile around sink i calculated 
using equation (pQ), dr indicates that the direction of integration is radial. Hence, we are essentially 
calculating the column density of gas along a line of length 2000 AU starting at the point of interest 
and pointing away from the sink. The sink particle that gives the highest column density at the 
point of interest is used to calculate the column density and also the local density at that point 
within the simulation. The local density is calculated using the density profile of the chosen sink. 
The choice of highest column density is made to ensure that we are not calculating the column 
density from the tail of the distribution around a more distant sink particle (see UED09). The 
limit of integration was set at 2000 AU to agree with the value used in UED09. 

As in UME10, we impose a minimum temperature T m ; n to the gas. Hence, if the calculation 
of T gas produces a value lower than T m i n , we use T gas = T m ; n instead. We consider the effects of 
changing the value of T m i n . 



3. The Simulations 

We performed two new simulations, G05 and G10, which include all the dust energetics dis- 
cussed in H2.1\ but also add a more accurate treatment of the gas energetics described in more 
detail in UED09. In the following results section, we will refer to these new simulations as using 
a "complete energetics" algorithm. For comparison, we also include two simulations, 105 and D05, 
that were presented in UME10. Our initial conditions are identical to those described in UME10. 
The parameters of the simulations are given in Table [TJ In the third column, Tj ust and T gas refer 
to the dust and gas temperature calculated using the complete energetics algorithm, while Tk is 
the actual temperature that we use for the gas. The letters I, D, and G stand for "isothermal," 
"dust," and "gas," respectively. 

In all simulations, the initial density of our cloud is p = 4.75 x 10~ 20 gcm -3 , or n = 1.22 x 
10 4 cm~ 3 assuming \x = 2.33. This density is simil ar to the media n average density (n = 1.6 x 10 4 ) of 



a well-studied sample of massive dense clumps (|Wu et alj|2010l ). In simulations 105, D05, and G05 



the minimum temperature of the gas is set at T m ; n = 5K, which corresponds to an initial Jeans mass 

Table 1. Numerical Parameters of the Simulations 



Run 


T min [K] 


T K 


Aftot [M e ] 


Lb ox [pc] 


M} nit [M ] 


Mj [M Q ] 


r acc [AU] 


105 


5 


T ■ 

± mm 


671.4 


0.984 


0.617 


0.0080 


152 


D05 


5 


^dust 


671.4 


0.984 


0.617 


0.0080 


152 


G05 


5 


T 

± gas 


671.4 


0.984 


0.617 


0.0080 


152 


G10 


10 


T 

± gas 


1898.0 


1.390 


1.744 


0.0226 


215 
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Mj mt = O.617M . Sink particles are created at a threshold density of p c = 2.822 x 10~ 16 gcm~ 3 , 
or n c = 7.252 x 10 7 cm -3 , which represents a contraction by a factor of 5942. The corresponding 
Jeans mass at p = p c is Mj = O.OO8M0. As in UME10, we adjust the particle mass such that 
each Jeans mass is resolved with 200 particles. The simulation starts with 64 3 particles, but we 
allow 2 levels of particle splitting (N gen = 2, see UME10). The effective number of particles is 
therefore 256 3 , the total mass of the system is M tot = (256 3 /200)Mj = 671.4M , and the box size 
is Lbox = (M to t/ p) 1 / 3 = 0.984 pc. Our SPH code uses a standard cubic spline smoothing kernel. 
The individual smoothing lengths are adjusted dynamically so that each gas particle has about 50 
neighbors. 

All simulations use identical initial conditions. Particles are laid down on a 64 3 cubic grid, and 
displaced in order to reprodu ce a Gaussian random density fluctuation with a P(k) oc k~ 2 power 



spectrum (IKlessen et al.lll998l . MES06). Initial velocities are then adjusted in order to reproduce a 
pure growing mode (MES06). In simulation G05, we set the minimum temperature to 5 K to allow 
a direct comparison with simulations 105 and D05 taken from UME10. However, most observations 
suggest Tk > 10 K, except in well-shielded areas of dense cores. Although massive, dense clumps 
may in fact be quite cold before they form stars, we explore the effects of T m \ n with a second 
simulation with full energetics, G10, in which we used T m ; n = 10 K. We did not change the initial 
density p and the threshold density p c . Doubling the temperature while keeping p c fixed increases 
the minimum Jeans mass by a factor of 2 3 / 2 , up to 0.0226M Q . Following the approach used in 
UME10, we decided to keep the same resolution for all simulations: 200 particles per Jeans mass. 
As a result, the particle mass and total mass increase by a factor or 2 3 / 2 , and the box size increases 
by a factor of 2 1 / 2 . The total mass Mtot, box size L DOX , the initial Jeans mass Mj mt , Jeans mass 
Mj at sink formation, and accretion radius r acc are listed in Table [H 

In all simulations, the first particle splitting occurs at a density p = 5.80/3, the second one at 
density p = 371 p, and sink formation starts at density p = p c = 5942p. The first sinks formed have 
an initial mass M ~ Mj. Sinks formed afterward will tend to be initially more massive since gas 
heating by the first sinks increases the Jeans mass. The initial free-fall time in our simulations is 
iff = 9.64 x 10 12 s = 3.06 x 10 5 yr. We run our simulations for a few free-fall times until the most 
massive sink particle in the simulation reaches M ~ 21M Q . We halt the simulations at this point 
since the luminosity from massive stars will produce signific ant ionizing phot ons. These photons 



will then dominate the evolution of the simulation as seen in lDale et al.1 (120051 ) . 



4. Results 

4.1. Fragmentation and Sink Formation 

Figure [T] shows the fractional mass distribution of the gas and sink particles as a function of 
time for the various simulations. The mass fraction in gas and in sinks is nearly identical for all 
runs with feedback (D05, G05, and G10). In all cases, the transition from gas- to sink-dominated 
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Fig. 1. — Evolution of the mass fraction in gas and sinks, and the number of sinks, for all four 
simulations. Top panel: run 105; second panel: run D05; third panel: run G05; bottom panel: run 
GIO. Top solid lines show the mass fraction in gas. Bottom solid lines show the mass fraction in 
sinks. Dotted lines with the scales on the right axes show the number of sinks. The scales of the 
right axes in the first and fourth panels are different from those in the second and third panel. The 
one in the fourth panel is larger than the one in the third panel by a factor of 2 3 / 2 to account for 
the increased volume. 
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mass fraction takes place around t = 2.45%, while for run 105 this transition takes place earlier, 
at t = 2.36%. By t = 2.4%, we have formed 74 sinks in run D05 and 118 in run G05. The 
slower formation of sink particles in run D05, for which Tk = 7d us t is due to the fact that the dust 
temperature is usually larger than the gas temperature, and setting the gas temperature to that 
higher value tends to inhibit the formation of sinks. We discuss this in more detail in 



It is interesting to compare runs G05 and G10. As we explained in §3, the volume simulated 
in run G10 is larger than the one simulated in run G05 by a factor of 2 3 / 2 . Hence, because run G05 
formed 118 sinks by t = 2.4%, we would expect run G10 to form about 334 sinks just because of 
its increased mass, not accounting for the effect of temperature. Run G10 actually forms 365 sinks 
by t = 2.4%, within 8% of our prediction. The right axis in the bottom panel of Figure [TJ has been 
rescaled by a factor of 2 3//2 compared to the third panel to allow a direct visual comparison. Sinks 
form faster in run G10 than run G05, but eventually the number of sinks per unit volume becomes 
comparable. Changing the minimum temperature affects the formation of the first sinks. But once 
several sinks have formed, the minimum temperature becomes irrelevant, because feedback heating 
by sinks raises the gas temperature above 10K in regions where the next sinks form, as we show in 
I below. 



In Table EJ we compare the final state of all four simulations. The third, fourth, and fifth 
column give the number of sinks, maximum sink mass, and mass fraction in sinks at the final time, 
respectively. Comparing runs D05 and G05, which both have a total mass of 671Mq, we find that 
fewer sinks are formed in run D05, but the gas that is prevented from forming new sinks ends up 
accreting onto the existing sinks, so the total mass in sinks is roughly the same at the end of the 
simulations. 

Figure [2] shows the distribution of sink and gas particles in our simulation box at t = 2.4% for 
run G05. We chose this time because this is when a substantial fraction of the total number of final 
sinks have formed and there is also a significant amount of gas remaining, as seen in Figure [TJ The 
morphology of the gas shown in this figure is similar to that seen in similar figures in UME10. There 
are fewer sinks formed compared to the isothermal simulation 105, but there are more compared to 
the simulation with only dust heating D 05. The dis tribution of sink particles near the intersection 



of filaments is similar to what is seen in iBatd (|2009l ). 



Table 2. Final State for all Simulations 



Run 


%nal [%] 


-^sinks, final 


M sink , max [M ] 


/sinks, final 


SFRg 


105 


2.5 


3429 


0.50 


60% 


0.24 


D05 


2.5 


74 


20.8 


50% 


0.20 


G05 


2.4 


118 


20.8 


46% 


0.19 


G10 


2.4 


365 


24.0 


47% 


0.20 



X 



Fig. 2. — XY position plot of sinks and gas particles at t = 2.4% for run G05. Black and blue dots 
indicate gas particles. Blue dots are gas particles which have undergone one particle splitting. Red 
dots are sinks. The box is 0.984 pc x 0.984 pc. 
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Fig. 3. — Solid histograms: separations between each sink at the time of its formation and the 
nearest already existing sink at that time. Dashed histograms: separations between each sink and 
its nearest neighbor at the end of the simulation. Top panel: run G05; bottom panel: run G10. 
Dotted lines indicate the value of 2r acc . 
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The sinks form in close proximity. To quantify this, starting with sink #2, every time a sink 
formed, we calculated the distance to the nearest existing sink. The resulting distributions are 
shown as solid histograms in Figure There is a "natural" separation Ar = 2r acc , which comes 
from the method used for forming sinks. The algorithm converts gaseous spheres of radius r acc 
into sinks. When a dense region fragments into several sinks, the minimum separation is of order 
2r acc , corresponding to spheres of gas that are in contact. Also, when a new sink forms near an 
existing one, the separation tends to be of order 2r acc : the gas density profiles near sinks ( §4.3|) 
show that the density decreases with increasing radius; hence the gas located the closest to the 
existing sink will be the first to reach the threshold density p c . The large separations seen in 
Figure [3] correspond to sinks that were the first to form in a gas clump that was located away from 
all other clumps where sinks were already present. Although there is a physical reason for sinks to 
form close to one another, the particular value of 2r acc has no physical meaning. The parameter 
r acc is adjustable, and so is the density threshold p c . Only the combination r^ cc p c has a physical 
meaning, since it determines the Jeans mass Mj. We could have used a larger value of r acc and a 
correspondingly smaller value of p c . Sinks would have formed farther apart, but gravity would have 
brought them together, and the end result would have been essentially the same. In very dense 
regions, separations between sinks can reach values as low as 10 — 20 AU. They are clearly much 
closer than they were when they formed. 

Although some sinks move closer together, others move apart. The final distribution of separa- 
tions is broad and skewed (dashed histograms in Fig. [3]), indeed bimodal. By the end of simulations 
G05 and G10, the mean separation is 4740 AU and 4940 AU, respectively, with final median sep- 
arations of only 700 AU and 556 AU, respectively. The med ians are much smaller than even the 



projected median separations in nearby clusters of 14800 AU (iGutermuth et al.ll2009l ). The median 
separations in both simulations increased by a factor of 3 or more between t/tg = 2.2 and 2.4, 
so further evolution may lead to larger mean separations. In addition, the overall d istributions, 



especially for G10, do resemble the observations (e.g. Fig. 2 of IGutermuth et al.ll2009i ) in having a 
tail toward large separations. Given that the confusion limit of the observations is about 3000 AU, 
they would not resolve the very close pairs. Indeed, the peak of the distribution below 100 AU 
might be interpreted as binaries and multiple systems. 



4.2. Evolution of Luminosity and Star Formation Efficiency 



Extragalactic observers employ a relation between total far- infrared luminos ity and star for- 
mation rate to estimate star formation rates in dusty starbursts (IKennicuttl 1 1 998l ) : 



SFR(M yr- 1 ) = 4.5 x 10- 44 L FIR (erg s^ 1 ) = 1.7 x 10- 10 L FIR (L Q ). 



(3) 



For convenience, we express this K98 relation as SFR/L = 1.7 x 10 ~ 4 M ffi Myr" 1 ^ " 1 or L/SFR = 
5.9 x 10 3 Lq Myr/M . A more recent calibration using the IMF of iKroupal (|2002l ) yields a slightly 
slower SFR, or higher L/SFR = 6.9 x lO 3 L Myr/M . 
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Wu et al.l (|2005l ) have used this relation for massive dense clumps in our galaxy, finding a similar 
relation between Lfir and the line luminosity of dense gas tracers like HCN for individual clumps as 
was foun d in starburst galaxies by iGao fc Solomonl (|2004j), as long as Lfir is above about 10 45 _Lq. 
However, iKrumholz &: Tanl (|2007l ) noted that the most massive stars, and hence the luminosity, 
take a considerable time to build up during cluster formation. Using an analytical prescription for 
star formation, they plotted L/SFR versus time, finding values 1-2 orders of magnitude below the 
extragalactic relation for times less than 1 Myr. 

In Figure HJ we show how the total luminosity L* ot , total sink mass M* ot , time-averaged star 
formation rate (SFR) = M* ot /i, and L* ot /(SFR) vary with time in our simulations. The build-up 
of stellar mass and luminosity is quite rapid, but the luminosity lags the mass. The calculations 
were stopped at 2.5tg or 0.75 Myr, before the peak of star formation (see Fig. [4]). Both the sink 
formation rate and the mass accretion rate onto sinks are still increasing at the end of the simulation 
(see Figured]). At that time, L* ot /(SFR) lies a factor of 10 below the K98 relation. When the total 
luminosity first exceeds 10 5 L Q , the value of L* ot /(SFR) lies about a factor of 20 below the K98 
relation. The values of L\ ot , M* ot , and (SFR) are larger for run G10 simply because the volume 
simulated is larger. To allow a direct comparison with run G05, we plot the results for run G10, 
scaled down by a factor of 2 3 / 2 (dotted line). The results are very similar to the ones for run G05. 



These results confirm the prediction by IKrumholz fe Tanl 1)20071 ) that the K98 relation will 
underestimate the SFR for an individual cluster at early times, but the discrepancies are somewhat 
less than they found. If star formation has proceeded for a few free-fall times, or about 1 Myr, 
the K98 relation becomes better for an individual clump. However, there must still be sufficient 
dust to convert most of the luminosity into far- infrared radiation. Once i onizin g radiation turns 
on, Lfir may become a better tracer of the SFR. IVutisalchavakul Evansl (|2012l ) found agreement 
to a factor of 2 between SFRs calculated from Lfir and radio continuum emission for a sample of 
massive dense clumps. 

The last column of Table [2] shows the star formation efficiency per free-fall time, SFRg, defined 
in IKrumholz fe McKed (|2005l ). It is essentially the ratio /sinks, final Aff- The values are similar for all 



runs, even though the number of sink particles is higher in the new simulations compared to the 
simulation with only dust heating energetics. The values in Table [2] assume that all mass entering 
the sink winds up in the star. As noted earlier, this almost certainly overestimates the stellar mass 
by factors of 2-3, suggesting values of SFRg- of 0.07 to 0.10. Star formation efficiencies in massive 
dense clumps are not well-constrained, but are certainly much lower. Using the numbers for the 
mean Lfxr/M V i T for massive dense clumps from lWu et al.l (|2010l ) and equation (|3j), along with the 
free-fall time at the mean density of these clumps (0.27 Myr), yields values of SFRg- around 0.006. 
If Lfir underestimates the SFR in those clumps by a factor of 20, the values would agree better 
with those in Table O but more likely the SFRg is lower than found in our simulation. Overly 
fast star formation is a common feature of simulations that do not include means to slow down the 
star formation process . Even recent simulations with radiative feedback, turbulence, and outflows 
(jKrumholz et al.ll2012l ) produce very similar values of SFRg- to those in Table [2j 




Fig. 4. — Total luminosity (top left), total stellar mass (top right), average star formation rate 
(bottom left), and luminosity per average star formation rate (bottom right) as a function of time 
for runs G05 and G10. All panels show evolution of a different variable as a function of time, in 
both absolute years (bottom axes) and free-fall time (top axes). The top solid lines, bottom solid 
lines, and dotted lines correspond to run G10, run G05, and run G10 scaled down by a factor of 
2 3//2 , respectively. Also plotted in the bottom right panel as a horizontal line labeled Kennicutt98 
is the relation given in equation ([3]). 
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4.3. Density Profiles around Sinks 



We calculate the density profile around each sink at each time step. The density profile will 
change around a sink as it moves within the cloud, gathering gas particles and/or as gas particles 
accrete irregularly onto it. In Sj2] we defined the density parameters, n Q and a, of the density 
profile, equation ([I]). Figure [5] shows the average values of the density profile parameters, a and 
n Q , for all the sinks that formed in run G05. We only included sinks that have accretion rates 
M > lO _8 M0yr _1 , because sinks with lower accretion rates have too few neighboring gas particles 
to allow for an accurate determination of the profile. We summarize the results in Tabled together 
with the results for run D05, taken from UME10. The values are similar for both runs. This suggests 
that the density profile surrounding the individual sinks is not strongly affected by the energetics 
algorithm used. 

Groups studying young star-forming cores have calculated density profiles from observations. 
Comparing our average density profile values to those obseryation ally derived values, we find excel- 
lent agreement with our values of a and n Q . IShirley et al.l (|2002l ) studied Class cores and found 
(a) = 1.63 ± 0.33 and a typical value of a = 1.8 ± 0.1 if they ignored two sources with aspherical 
emission contours. I Young et al.l (|2003l ) studied Class I cores and found (a) = 1.6 ± 0.4. Our sim- 



ulations are consistent with either of these v alues for a. The va lues of log n derived from these 
two studies are (log(n /cm- 3 )) = 6.1 ± 0.2 (IShirlev et al.1 12002 1 and (log(ra /cm~ 3 )) = 5.4 ± 0.5 
( Young et al.ll2003l ). Our value of (log(n /cm~ 3 )) is 6.4 ± 0.3, including only points with accretion 
rates greater than 1Q~ 8 M^ yr -1 . These values tend to agree better with the results for the Class 
core study of I Shirley et al.l (|2002l ). suggesting that most of the sinks are reflecting early stages in 
star formation. 



4.4. Evolution of Temperature, Density and Velocity 

Figure [6] shows the density and temperature evolution of gas particles in our simulations. The 
behavior of the gas particles in these figures resembles what was seen in similar figures in UME10. 
For example, at early times the temperatures of the gas particles are confined to low values at all 
densities. However, as the simulation evolves and more sink particles form, the gas particles are 
heated and tendrils of gas particles appear in Figure [6j They are reaching toward high-density and 



Table 3. Density Profile around Sinks 



Run 


(a) 


(log(n /cm 3 )) 


D05 


1.7 ±0.4 


6.5 ±0.3 


G05 


1.7 ±0.4 


6.4 ±0.3 



-15- 



1.6 t„ 1.7 t„ 1.8 t ff 2.0 t„ 2.4 t„ 




20 40 60 80 100 120 

Sink Number 



Fig. 5. — Density profile parameters for all 118 sinks formed in simulation G05. Values of density 
profile parameters, a and n , are shown for the sinks. Error bars shown with a solid line indicate 
the standard deviation for each individual sink. Error bars shown with a dotted line indicate the 
minimum and maximum value that a and n Q take during the course of the simulation. 



-16- 




2468 10 2468 10 2468 10 
log(n) [cm -3 ] log(n) [cm -3 ] log(n) [cm -3 ] 




2468 10 2468 10 2468 10 
log(n) [cm -3 ] log(n) [cm -3 ] log(n) [cm -3 ] 



Fig. 6. — Temperature and density of gas particles as a function of time for runs G05 (top 6 panels) 
and G10 (bot tom 6 panels) , at various times, as indicated. The green line shows the equation of 
state given by [Larson] (120051 ). The horizontal and vertical red lines show the minimum temperature 
Tm in and threshold density n c for sink formation, respectively. The statistics of the temperature 
are listed in Table HI 
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high-temperature regimes. These tendrils show the behavior of the gas density and temperature 
near luminous sink particles. Another interesting feature seen in this figure is the slight increase 
of the floor of the temperature distribution as a function of time. This feature of an increasing 
temperature floor was also seen in run D05 (UME10). In Figure [6j we also plot the equation of 
state given by [Larson] (|2005l ). As in run D05 (UM E10), the gas does not follow this equation. This 
is not surprising because the equation of state of iLarsonl (120051) does not take into a c count stellar 



heati ng. This has been confirmed to be an important effect in lKrumholz et al.1 (|2007l ). lQffner et al 
(|2009h . and UME10. 



The horizontal and vertical red lines in Figure [6] show the minimum temperature T m i n and 
the density threshold p c , respectively. Gas particles located on the right of the vertical line, in 
regions where sinks form, all have temperatures significantly larger than T m [ n . These particles do 
not turn into sinks because their high temperatures make them unbound [see eq. (13) in MES06]. 
This explains the results that were presented in Figure [1] and Table [2j Increasing the minimum 
temperature from 5K to 10 K might affect the formation of the first few sinks, but once these sinks 
have formed and started reheating the gas, the formation of subsequent sinks is unaffected by the 
particular value of T m \ n , because the gas temperature in regions of sink formation is already much 
larger. This explains why the gas fraction drops at essentially the same rate in simulations G05 and 
G10 (Fig. [1]), and why both simulations form roughly the same number of sinks per unit volume. 
Judging from Figure [H it would take a minimum temperature of order 50 K or more to make a 
difference. 

Most of the cores in our simulations are located inside dense filaments. The few cores located 
in low-density regions actually formed inside filaments, and were later ejected by 3-body encounters 
(MES06). This concentration of cores can explain the large gas temperatures found in these regions 
(even though more cores means more competition for accreting gas from the same region, possibly 
leading to lower accretion luminosities). T he spatial dis t ribut ion of cores is a consequence of 
the assumed periodic boundary conditions. iHeitsch et al.l (|2008l ) performed simulations of cloud 
fragmentation with isolated boundary conditions, and showed that these initial conditions lead to 
stronger initial fragmentation in the early stages, and a more distributed core formation at late 
stages. Observations are, however, finding st rong concentrations of cores and protostars along 
narrow filaments (e.g.. iHennemann et al.ll2012l ). very similar to patterns seen in the simulations. 



TableU]shows various statistics for runs with heating (D05, G05, and G10), i.e., the average gas 
and dust temperatures, and average density (calculated by averaging over all gas particles), number 
of sink particles, percentage of mass in sink particles, and mass of the most massive sink particle, 
at different times. We find that the average gas and dust temperatures both increase throughout 
the simulations. This increase can be explained by the increasing number of sink particles and the 
increasing mass of the individual sink particles, both of which lead to more substantial radiation 
fields. 



Since it is essentially the mass of sinks and accretion rate of gas onto sinks that determine 
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Table 4. Particle Statistics for Runs with Heating 



Run 


time [iff] 


Tgas [K] 


T dust [K] 


log n [cm 3 ] 


-^sinks 


/sinks 


max(M sink ) 


D05 


1.6 


5.40 ± 1.91 


5.40 ± 1.91 


4.65 ± 0.67 


3 


0.02% 


0.07 


D05 


1.7 


9.28 ± 3.84 


9.28 ± 3.84 


5.02 ±1.08 


16 


0.4% 


0.72 


D05 


1.8 


14.8 ± 6.08 


14.8 ± 6.08 


5.12 ± 1.12 


28 


2.4% 


2.11 


D05 


2.0 


24.5 ± 9.66 


24.5 ± 9.66 


5.05 ± 1.15 


53 


12% 


7.99 


D05 


2.2 


33.0 ± 14.3 


33.0 ± 14.3 


4.88 ± 1.16 


70 


28% 


12.40 


D05 


2.4 


48.3 ± 29.7 


48.3 ± 29.7 


4.76 ±1.21 


71 


45% 


17.99 


G05 


1.6 


5.64 ± 1.84 


4.50 ±2.20 


4.63 ±0.68 


3 


0.03% 


0.06 


G05 


1.7 


7.61 ±4.14 


9.64 ± 3.89 


5.16 ± 1.15 


26 


0.5% 


0.78 


G05 


1.8 


10.8 ±7.22 


14.9 ± 5.94 


5.41 ± 1.20 


49 


2.8% 


1.7 


G05 


2.0 


20.6 ± 10.8 


25.9 ± 10.6 


5.40 ± 1.19 


89 


14% 


7.9 


G05 


2.2 


28.2 ± 18.1 


33.6 ± 17.8 


5.19 ±1.21 


108 


30% 


12.8 


G05 


2.4 


35.6 ±25.6 


43.0 ± 25.4 


4.89 ± 1.18 


118 


46% 


20.8 


G10 


1.6 


10.2 ± 1.67 


5.35 ±2.60 


4.65 ±0.67 


5 


0.03% 


0.12 


G10 


1.7 


11.4 ±4.24 


11.8 ±4.78 


5.21 ± 1.11 


37 


0.6% 


2.22 


G10 


1.8 


17.5 ± 7.07 


19.6 ±7.11 


5.63 ± 1.26 


138 


3.1% 


6.08 


G10 


2.0 


26.7 ± 12.8 


30.0 ±13.1 


5.76 ± 1.22 


301 


15% 


15.1 


G10 


2.2 


33.1 ± 18.2 


35.7 ± 18.0 


5.57 ± 1.22 


346 


32% 


21.8 


G10 


2.4 


39.6 ±22.6 


42.7 ±22.5 


5.41 ± 1.24 


365 


47% 


24.0 
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the dust temperature, we expect the temperature of the gas being accreted to have little effect if 
it accretes supersonically, which is indeed the case (see Table [5] below). At all times t > 1.7tg 
in both G05 and G10, the mean gas temperature is slightly cooler than the dust temperature. 
This is expected since some gas is at quite low densities (n ~ 10 3 or 10 4 cm -3 ) and the rest is no 
hotter than the dust temperature. The gas is fully equilibrated to the dust temperature only at 
quite high densities, such as n > 10 6 cm~ 3 . The dust temperature is initially fixed at 5K for run 
G05 and 10 K for run G10, and it can change only when heating via nearby stellar radiation fields 
begins. Figure [7] illustrates these effects in more detail. It shows the histograms of the dust and gas 
temperature at various times for low-density (n < 10 5 cm -3 ) and high-density (n > 10 5 cm -3 ) gas. 
Except for the first time step shown, the dust temperature is higher than the gas temperature at 
low densities (red histograms). It also shows that the dust and gas temperatures are nearly equal 
at high densities (black histograms) due to dust-gas collisional coupling. At t = 1.6%, the dust is 
cooler on average than the gas, for low density gas (see also, the first line of Table 0]). At these 
early times, the dust has not yet been heated to high temperatures by forming stars. Because the 
dust temperatures are lower, the dominant heating source for low-density gas in our simulations is 
cosmic-rays, which can raise the gas temperature above the dust temperature. 

Figure [8] shows the dust-gas temperature difference as a function of density at various times 
for run G05. The contours indicate the region in parameter space, Td ust — T gas versus log n, where 
the majority of the material is located. The gas is confined to the same region for most of the 
simulation, except during early times when the dust is heating up due to formation of luminosity 
sources. As seen before in Table [4] and Figure [7] during the earliest times, the gas temperature 
is hotter than the dust temperature because of the low-intensity of the radiation field. At high 
densities (log n > 8), the SPH gas particles tend to approach the horizontal line which indicates 
equal dust and gas temperature. At these high densities when the dust and gas temperatures are 
comparable, dust heating dominates all other forms of heating/cooling for the gas. This result was 
seen in UED09. By 2tg, almost all the gas has Tk within 30% of Td ust , and a large fraction is even 
closer to equilibration. Figure [8] also shows that high dust temperatures are most likely to exist in 
high-density regions (as seen by the points), presumably due to the proximity of a nearby forming 
star. 

Figures [7] and [8] show that neither the gas nor the dust are well described by isothermal or 
barotropic equations of state once protostars begin to heat their surroundings. There is a strong 
effect on the thermal behavior of the clump once protostars develop significant luminosity. That 
effect occurs during the sink formation process. 

For each gas particle accreted onto a sink, we calculated the velocity v of the particle relative 
to the sink, and the Mach number A4 = v/c s , where c s is the sound speed at the location of the 
sink. Table [5] shows the mean values of v and M, and also the mean values of the accretion rate 
M, for all simulations with heating. Again, similar behavior of the sinks is seen for both G05 and 
G10, and also D05. On average, particles are being accreted supersonically. The sound speeds near 
sinks are comparable for runs G05 and G10, since, as we saw, the temperature in these regions 
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Fig. 7. — Dust and gas temperature histograms in the low and high density regimes, for run G05. 
The dust (solid line) and gas (dotted line) temperature histograms are plotted for low (red) and 
high (black) densities. 
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Fig. 8. — Temperature difference versus density shown at various times for run G05. Horizontal 
line marks equal dust and gas temperature. Contours show 3, 10, and 20 a contour levels of density 
of SPH gas particles. Dots indicate gas particles which have dust temperatures greater than 100 K. 
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greatly exceeds T m i n . Hence, the slightly larger accretion velocities for run G10 results in slightly 
larger Mach numbers. Compared to the runs with full energetics, run D05 has on average larger 
accretion rates, larger accretion velocities, and comparable Mach numbers. This run forms fewer 
sinks, but the total mass in sinks is the same (see Fig. Q] and Table [2]). Hence, the same mass 
is being accreted by fewer sinks, resulting in higher accretion rates. Also, fewer sinks means that 
they are on average more massive. Since the accretion velocities are highly supersonic, gas particles 
accrete essentially at the escape velocity, which is larger for more massive sinks. However, the gas 
temperature, and consequently the sound speed, is also larger for run D05. As a result, the Mach 
numbers are comparable to the ones for runs G05 and G10. 

These infall velocities are considerably higher than the assumed turbulent broadening (Doppler 
b parameter is set to 1 km s _1 as discussed in §2.ip . To the extent that these larger velocities allow 
for greater escape probabilities for photons in the primary cooling lines, the gas could be somewhat 
cooler than we calculate in the infalling gas around a sink. On the other hand, these are just 
the regions where the density becomes large enough to overwhelm gas cooling and couple the gas 
temperature to the dust temperature, so the effects are probably minimal. 



4.5. Mass Functions 



T he Ga l actic F ield Star IMF has been studied by many groups (see lSalpete 3l955l : lMiIler fc Scald 



1979: 


Scalo 


1986: 


Kroupa 


2002 





Chabrierl 120031 ). The cluster IMF has also been studied in 



young clustered, star-forming regions such as Ori on (Hillenbrand 19 97J) and in m ore isolated star- 
forming regions in t he Taurus (IBriceho et al.ll2002l ). Lupus (jComeron et al.ll2009r), and C hamaeleon 
( Alcala et al. 1997 ) molecular clouds. In a recent review paper, Bastian et al. (hold ) found no 
significant variations in the IMF for present-day star formation. 

In large-scale simulations of cluster formation, various treatmen ts of thermal energ etics were 
used. The earliest attempts assumed an isothermal equation of state (jKlessen et al.lll998l : MES06). 
In the work of MES06, we found that a properly-resolved isothermal calculation could only produce 
very low-mass fragments. The mass function in these simulations was log-normal with a peak 
that was determined by the resolution limit of the simulation. Other simulations attempted to 



Table 5. Statistics of Accretion onto Sinks. 



Run 


M 


v [kms 1 ] 


M [1O~ 5 M yr" 1 ] 


D05 


10.8 ± 7.4 


9.2 ±8.2 


2.20 ±3.08 


G05 


8.2 ±5.9 


5.5 ±5.3 


1.65 ±2.32 


G10 


10.8 ±7.6 


6.5 ±6.0 


1.37 ± 2.36 
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address this issue by modifying the equation of state to account for the increase in optical depth 
in higher density gas. They used a barotropic equation of state i n which the temperature of the 
gas increases at higher densities (IBate et al.l 1 2003; Li et al.l 120031 ; iBatd 120051 : IJappsen et al.ll2005l ; 
Larson! 120051 ; iBonnell et al.l 120061 ; IClark et al.l l2008l ) . In more recent work, groups have included 
various treatments of ra diative transfer to mo r e accurately account for the heating effect of the 



young stars on the gas (jKrumholz et al 



2010 



Bate 



2009 



Qffner et al 



2009 



Smith et al 



2009 



UME10). Inclusion of the heating effect shifts the mass function to higher masses. 

Figure [9] shows the mass functions of the sink particles at various times in the simulations G05 
and G10. At early times, the mass function is dominated by accreting low-mass objects. As the 
simulations evolve, these low-mass sinks continue to accrete and become more massive. At later 
times, the formation of low-mass sinks become less frequent because the gas is now hotter and the 
Jeans mass is higher, causing a delay in sink formation until they have reached the new higher 
Jeans mass. Although there is a tendency for the mass function to shift its peak to higher masses 
with time, the IMF in simulation G10 becomes relatively stable after about 1.8 iff. The shift to 
higher masses relative to the isothermal calculation agrees with the results of the references given 
above. 



We also plot the IMF's from ISalpeterl (|1955luGhabrieri (|2003h . and iKroupal (|2002h in the last 
panel for each run. We normalize the mass functions using the maximum mass in the simulation 
(solid curves) and the total number of objects (dashed curves), both given in Tabled The minimum 
mass used to create the IMF's was O.O8M and O.O1M for the ISalpeterl (|1955l ) and IChabrier 



(|2003l ) / iKroupal (J2002J) IMF's, respectively. Normalizing the mass functions based on the total 
number of objects highlights the fact that we form too few objects given the mass of the most 
massive sink in each simulation. The IMF's normalized to the maximum mass demonstrate that 
we can roughly reproduce the slope of high-mass tail above masses of ~ 1M & . 

As in UME10, our mass functions slightly over-predict the number of high-mass objects and 
miss a substantial fraction of low-mass objects. In UME10, we predicted that by including molecular 
cooling and cosmic-ray heating we would see a decrease in high-mass objects and an increase in 
intermediate- to low-mass objects. This has turned out to be the case and can be seen in the top 
panel of Figure [10] where we show the mass functions for our three different simulations with a 
minimum temperature of 5K. We form more intermediate-mass objects, while still retaining the 
slope at higher masses. 

In Table El we give the average and median sink masses in our simulations. We also give these 
values predicted by various empirical IMF's. We use a lower mass limit of 0.08Mq for the Salpeter 
IMF and O.OIMq for the Kroupa and Chabrier IMF's. Including the complete energetics algorithm 
(G05 and G10) decreases the resulting mean mass, compared to D05, as predicted. The change in 
the median mass shows a similar behavior. Although our median and average masses decrease, they 
still give values which are in disagreement with the empirical IMFs. Next, we perform a thought 
experiment to address this issue. 
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Fig. 9. — Mass histograms for runs G05 (top panels) and G10 (bottom panels), shown at different 
times. The tim e and number of sinks are indic ated in each panel. The analytic mass functions of 



Salpeterl (j 19551 ) (straight line) , IChabrierl (120031 ) (curved line) , and iKroupal (|2002l ) (segmented line) 



are shown as solid lines when normalized to maximum mass object, and dashed lines when normal- 
ized to total number of objects. Vertical dashed lines show initial Jeans mass in the simulation. 
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Fig. 10. — Top panel: mass histog rams for all simulations with Tmm — 

5K. Mass function are 

shown at f inal time (2.5fo r for 105 and D05, 2.4% for G05). The empirical IMF from IChabrier 
(|2003l ) and iKroupal (|2002l ) are shown, normalized to a maximum mass of ~ 20M Q . Bottom panel: 
mass histograms for runs D05 and G05, and efficiency of 30%. 
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The regions around our sinks have density profiles similar to those of gas around Class objects 
discussed in £j4.3l Thus, it is reasonable to associate them with dense cores. As this material falls 



into the sinks, we have so far assumed that all the mass winds up in a star. lAlves et al.l (|2007l ) claim 
that the dense core mass function is related to the stellar IMF, i. e., sharing the sam e shape, but 



shifted to higher masses. This is interpreted as a core-to-star efficiency factor of 30%. I Enoch et al 



(|2008l ) find a lower-limit of 25% for the core-to-star formation ef ficiency. The missing mass is likely 



to be removed in molecular outflows driven by stellar winds (e.g.. lDunham et al. 



2010; 



Hansen et al 



2012 ) . These outflows are not included in our simulation. To include this inefficiency, we scale the 



mass of sinks in the data from our simulations down by a factor of 0.30. We then recalculate 
the shape of the empirical stellar IMF for the new maximum masses and total number of stars. 
We show the results in the bottom panel of Figure [TU1 We plot the empirical IMF normalized 
to the new maximum mass. Our results exhibit better agreement with a stellar mass function 
with a 30% efficiency factor. Since we have only decreased the masses, the high-mass slope in 
our new simulation data remains the same and agrees with the empirical IMF. In addition, the 
peak in the IMF shows better agreement. At 0.1M Q , our new mass function (for the complete 
energetics algorithm) is only missing a factor of 2-3 objects compared to the empirical IMF, which 
is significantly less than the factor of ~ 20 seen in Figure [TO} However, the adjusted mean and 
median masses are still higher than in the empirical IMFs (Tabled]). Since we assumed that all 
the mass of the sink went into the star in order to calculate luminosities, it is not self-consistent 
to apply the efficiency factor after the calculation. We do it only to suggest avenues for future 
exploration. There are qu estions about the connection betwee n the core mass function and the 
stellar mass function (e.g., IClark et al.l 120071 ; ISmith et al.l 12009 ). Other effects may decrease the 
mass of the final stars, including di sk fragmentation to form binary or multiple stars and brown 
dwarfs (e.g.. IStamatellos et al.ll201ll ). 



5. Summary and Conclusion 

The motivation for this work was to explore the effect of dust-gas energetics in a clustered 
star formation simulation. We have presented the results of two new simulations (G05 and G10), 
which include the complete energetics algorithm discussed in UED09. The properties of the models, 
including two previous ones (105 and D05) are summarized in Table [H and we compare them here. 

As in D05 and 105, sink particles, representing protostars, form along filaments and especially 
at intersections of filaments. The temperature of a simulation which included gas cooling was on 
average lower than in a simulation with Tk = 7d us t, as expected. With a lower average temperature, 
fragmentation was more prevalent and more objects were able to form. The average density profile 
parameters surrounding a sink were similar among the four simulations and agreed with observations 
of low-mass Class sources. However the infall speeds were significantly supersonic, and mass 
accretion rates were high, both in contrast to observations of low-mass protostars. Infall speeds 
and mass accretion rates for high-mass protostars and protostars in clustered environments are 
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poorly constrained. Infall speeds and mass accretion rates are somewhat smaller for G05 and G10 
compared to D05. 

We added a calculation of the ratio of far- infrared luminosity (£fir) over t ne SFR to test 
the use of .Lfir as a SFR tracer in very young regions of clustered star formation. We found 
that Lfir/SFR increases rapidly during the simulation, but that it is significantly lower (factor 
of 10) than the ratio used to measure extragalactic SFR at the end of the simulations (around 
0.7 Myr). Measurements of SFR for very young clusters (ages < 1 Myr) using Lfir are very likely 
underestimated. 

We computed the mass evolution of protostars during the simulation, and we compared the 
mass function at the end of each simulation. In our previous work, UME10, we found that a non- 
isothermal, stellar-source dependent energetics algorithm radically reduced the number of young 
stars that were formed and formed more massive stars, compared to simulations with isothermal 
gas. However, the simulations in UME10 over-produced high-mass objects and missed a large 
fraction of low-mass objects. We predicted in UME10 that including a more realistic calculation of 
the gas temperature might address this problem. In this work, we included the complete dust-gas 
energetics algorithm. This change increased the number of intermediate mass objects, but the 
deficiency of low-mass objects persists. 

The two main differences between D05 versus G05 and G10 were the temperature distribution 
and the mass function, which are related to each other. In a lower temperature environment 
with more sink particles forming, there was less material available to be accreted and therefore 
a smaller fraction of massive objects were formed. This affected the mass function and led to a 
slight decrease in the number of high-mass objects and an increase in the amount of low-mass 
objects when compared to the simulations with dust heating only in UME10. We found very little 
difference in the mass function between G05 and G10, indicating that the initial temperature is not 
very important; feedback from the first protostars rapidly erases the effects of initial temperature. 

We performed a thought experiment in which we tried to explain the discrepancy between our 
mass function and the empirical IMF. In our simulation, we assumed that all the mass in a sink 



particle winds up in a single star. However, studies of nearby clouds (e.g., lAlves et al.l 120071 and 



Enoch et al.ll2008l ) show that about 70% of the core mass does not wind up in the star, probably 



becau se it is removed by stellar winds and the resulting molecular outflow (e.g., iDunham et al 



2010 ) . If we multiply our sink masses by 0.3, we get better agreement with the IMF. Although we 
have shown that including dust-gas energetics is essential, other effects (e.g., magnetic fields, turbu- 
lence, outflow s , etc.) will need to be included for a full understanding. In a promising development, 



Hansen et al.1 (|2012l ) found that including outflows along with radiative feedback reduced the mass 



accretion rates and protostellar mass es, hence luminos i ties, a llowing more fragmentation and better 



reproducing the IMF. More recently, iKrumholz et al.l (|2012l ) have found better agreement with the 



IMF when turbulence, outflows, and radiative transport are included, although their IMFs are still 
somewhat top-heavy. 
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Table 6. Average and Median Masses 



Run / observations 


Median 


Average 


D05 


2.26 


4.63 


G05 


1.63 


2.65 


G10 


1.50 


2.46 


G10 (x0.30) 


0.45 


0.74 


Sabeter C1955) 1 


0.13 


0.26 


Chabrier (2003^ 


0.10 


0.31 


KrouDa ('2002) 1 


0.15 


0.34 



Median calculated using normalization 
to total number of objects formed; Average 
calculated using normalization to maximum 
mass object. 



